arXivrl 509.03349V 1 [q-bio.MN] 10 Sep 2015 


Control of Stochastic and Induced Switching in Biophysical Networks 

Daniel K. Wells,William L. Kath,i'2.3 ^nd Adilson E. Motter^’S^.g 

^Department of Engineering Seienees and Applied Mathematies, 

Northwestern University, Evanston, IL 60208, USA 
‘^Northwestern Physieal Seienees-Oneology Center, 

Northwestern University, Evanston, IL 60208, USA 
^Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA 

^Department of Physies and Astronomy, Northwestern University, Evanston IL, 60208, USA 

Noise caused by fluctuations at the molecular level is a fundamental part of intracellular processes. 
While the response of biological systems to noise has been studied extensively, there has been limited 
understanding of how to exploit it to induce a desired cell state. Here we present a scalable, quanti¬ 
tative method based on the Freidlin-Wentzell action to predict and control noise-induced switching 
between different states in genetic networks that, conveniently, can also control transitions between 
stable states in the absence of noise. We apply this methodology to models of cell differentiation 
and show how predicted manipulations of tunable factors can induce lineage changes, and further 
utilize it to identify new candidate strategies for cancer therapy in a cell death pathway model. 
This framework offers a systems approach to identifying the key factors for rationally manipulating 
biophysical dynamics, and should also find use in controlling other classes of noisy complex networks. 

Subject Areas: Complex Systems, Biological Physics, Nonlinear Dynamics 


I. INTRODUCTION 

Cellular systems are not entirely deterministic, but 
are instead impacted by small, random fluctuations in 
the number and activity of molecules of intracellular 
species Such fluctuations lead to macroscopic ef¬ 

fects in a diverse array of processes. In differentiation, 
the resulting noise plays a central role in cell fate deter¬ 
mination and can allow clonal populations of differenti¬ 
ating cells to achieve distinct final states Hi. Noise can 
also produce spontaneous transitions, whereby it causes 
a system to switch from one stable state to another, often 
producing a significant change of phenotype or function. 
Such stochastic state switching occurs, for example, in 
the lac system, where rare, brief transcription events in 
the “off” state cause large bursts in LacY expression, 
which in turn can be amplified and stabilized by a posi¬ 
tive feedback loop [5[ . Stochastically-induced transitions 
also underlie recent observations of spontaneous dediffer¬ 
entiation in cancer cells [H, 0 , in which cancer stem cells 
arose de-novo from non-stem cell populations. 

The response to noise and the overall behavior of many 
biophysical systems are determined by an underlying epi¬ 
genetic landscape [8[. In this landscape, the valleys rep¬ 
resent the distinct achievable states of the system and 
the heights of the separating barriers determine their ro¬ 
bustness to noise. A benefit of this representation is that 
bifurcation points—locations in the parameter space at 
which one or more stable states suddenly cease to exist— 
correspond precisely to the points where one or more of 
such barriers first reach zero height as parameters change. 
This landscape thus incorporates two distinct features of 
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a state, namely its robustness to noise and its determin¬ 
istic stability, into one: the less robust a state is to noise, 
the closer it is to being eliminated through a bifurcation, 
and vice-versa. 

The landscape representation has been given a quanti¬ 
tative foundation as the quasipotential of the determin¬ 
istic component of the system dynamics [9[ and has been 
explored in experiments—e.g., to show how two parame¬ 
ters in the yeast galactose signaling network, the con¬ 
centrations of galactose and intracellular GalSOp, can 
alter the rates of stochastic switching in this bistable 
circuit 0- Despite these advances, researchers’ abil¬ 
ity to control this landscape in order to induce prespeci¬ 
fied biological outcomes has been generally limited to at 
most two parameters [ullil, and no general method 
currently exists to systematically tune transitions be¬ 
tween stable states and/or eliminate undesired states al¬ 
together. The possibility of such control would offer clear 
opportunities. For example, under the widely supported 
stochastic model for induced pluripotent stem cell gen¬ 
eration [TsI, a majority of cells have the possibility of 
being reprogrammed, even though existing technologies 
have achieved substantially smaller yields . The abil¬ 
ity to control the response to noise of differentiated and 
stem-cell states (e.g., inhibiting transitions to the first 
and promoting transitions to the second) could lead to 
enhanced procedures to create induced pluripotent stem 
cells. Similarly, in the context of the “cancer attractor” 
hypothesis [l5|, , in which normal and cancer cells cor¬ 

respond to distinct co-existing stable states, identifying 
interventions that destabilize, or eliminate, the cancerous 
state could lead to new therapeutic strategies. 

In this paper we propose a broadly applicable method, 
here termed optimal least action control (OLAC), that 
can predict and control the dynamical behavior and re¬ 
sponse to noise in a wide class of biophysical networks. 
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FIG. 1. Control of the response to noise illustrated for a multi-well quasipotential. (a) Without noise the state of the system 
is fixed and will not change over time, (b) In the presence of noise the state of the system wanders within the attraction 
basin and will be eventually ejected into the attraction basin of another stable state. For small noise such transitions are 
exponentially more likely to occur through the lowest barrier, even though other transitions are also possible in principle, (c) 
Using OLAC on a set of tunable parameters, we can alter the quasipotential to produce a desired response to noise, in this case 
tailored to increase transitions from state 2 to state 1 and to reduce transitions to state 3. (d) If desired, OLAC can also find 
combinations of tunable parameters that can alter the quasipotential in order to eliminate a stable state through a bifurcation; 
in this illustration state 2 is eliminated, in favor of state 1. (e-h) NEST for each of the quasipotentials in (a-d), respectively, 
where the nodes represent stable states and the continuous edges represent transition rates. A wider continuous edge indicates 
a higher transition rate. The dotted edges in (e) indicate transitions that could occur in the presence of noise. 


As schematically illustrated in Fig. [TJ the essence of our 
approach is that to control a biophysical system it is suffi¬ 
cient to identify interventions—e.g., changes to gene ex¬ 
pression, protein levels, or interaction rates—that can 
reshape the topography of the underlying quasipoten¬ 
tial in a desired way. This approach ultimately leads 
to a network of state transitions (NEST) describing the 
transitions between stable states and that can be con¬ 
trolled by changing the heights of the separating barriers 
without changing any quality of the noise. For a given 
system, this is achieved by determining the minimum 
action paths—those followed by the most likely noise- 
induced transition trajectories—and the corresponding 
transition rates between all pairs of stable states, and 
then optimizing these transition rates for a desired out¬ 
come. Furthermore, this general foundation in a physical 
least-action principle allows OLAC to be applied broadly 
to many other complex networks as well. In particular, 
while we focus our application of OLAC to biophysical 
networks, applications to other networks where noise and 
multistab ility play important roles, including power-grid 
networks [I 3 , polymer networks [l8| and food-web net¬ 
works [ 1 ^, among others, are immediate within the for¬ 
mulation we establish here. 

We apply OLAC to several gene network models and 
illustrate how this method can be used to make biologi¬ 
cally realizable reprogramming predictions. In the limit 
of zero noise intensity, OLAC automatically identifies bi¬ 
furcations that eliminate undesirable states and induce 
purely deterministic transitions to the desired ones. The 
significance of the latter is demonstrated by considering 
a third application, to eliminate cancerous states in a 
cell death network model, which concerns a time scale 


for which stochastic switches can be neglected. As il¬ 
lustrated in these examples, the NEST is a powerful yet 
simple representation that captures the essence of the 
state switching dynamics and can inform counterintu¬ 
itive results—e.g., the possibility of transitions through 
intermediate stable states when direct transitions are es¬ 
sentially impossible (a behavior observed even in high 
dimensions, where indirect transitions generally require 
longer paths). The method proposed here is easily im- 
plemetable and the computational effort scales linearly 
in the number of control parameters and the dimension 
of the state space, allowing our approach to be applied to 
large networks and high-dimensional systems in general. 


II. THEORY 

A. Transition Rates for Small Noise 

We consider biophysical networks whose determinis¬ 
tic components are described by non-linear differential 
equations of the form dXjdt = F(X;f2), where X is a 
vector representing the activity of the relevant biologi¬ 
cal factors, F is the function representing the rates of 
change of these factors, and is a set of tunable pa¬ 
rameters, which we show can be manipulated to drive 
cellular processes in advantageous directions. We fo¬ 
cus on the most prevalent case of systems with two or 
more stable states and, although our approach is gen¬ 
eral, for concreteness we first assume that these states 
are time-independent. Time-independent stable states 
correspond to fixed points X* = X*(f2), defined by 
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= 0, towards which neighboring trajectories 
converge over time. The set of all stable states represent 
the possible long-term behaviors of the deterministic sys¬ 
tem. 

Stochasticity is modeled here as additive Gaussian 
white noise, 

dX = F{X; n)dt -h A/i dW, (1) 

where 5 is the variance of the distribution (^her cases are 
discussed in the Supplemental Material [2^). With the 
addition of this small noise term, trajectories no longer 
approach the stable states asymptotically as in the deter¬ 
ministic case. Instead, a trajectory close to a stable state 
will oscillate stochastically within its basin of attraction, 
typically staying close to the fixed point for long peri¬ 
ods of time. The trajectory will also make rare but large 
excursions from the stable state. After sufficiently long 
time an excursion large enough to eject the trajectory 
from the original basin of attraction will necessarily oc¬ 
cur, at which point it will transition to the neighborhood 
of another stable state of the system. The time scales 
for the occurrence of such transitions may be shorter or 
longer than the biologically relevant ones. The manipu¬ 
lation of these time scales underlies much of the control 
approach introduced below. 

For a given noise intensity 5, the transitions between 
two stable states, i and j, occur as a Poisson process with 
a certain rate Rf j{Cl). These rates can be computed by 
evolving Eq. o, but in general at an unreasonably high 
computational cost. As a way to reduce this effort, we 
employ an asymptotic formula HH: 

Rlj (^) oc exp > (2) 

where ■ serves as an excellent approximation to Rf j 
for 5 small compared to 5'*^, which is typically the case 
for noise associated with biophysical systems. Here, S^j 
is the minimum of the Freidlin-Wentzell action 5'['], a 
functional over all possible transition paths (^i^j connect¬ 
ing the two stable states in the state space. The path (j)* j 
minimizing the action for the given pair of stable states is 
calculated numerically employing an implementation of 
the adaptive minimum action method [^, in which we 
determine this minimum action path through the opti¬ 
mization of a discretized version of the action functional 
using a quasi-Newton method (Appendix!^. Conceptu¬ 
ally, S^j represents the cumulative height of all saddle 
points traversed by the minimum action path between 
the states i and j. It should be noted that a propor¬ 
tionality constant is omitted in the expression for Rf j. 
Although there is no known formula for computing this 
constant in general [2^, the key condition for neglecting 
it is clear: as long as the actions Sij associated with two 
paths differ by more than 0{e)^ the exponential term 
will dominate and the omission of the proportionality 
constant will not affect the result qualitatively. 


The rate Rf j represents the transition probability per 
unit of time along the most likely direct path(s) between 
the stable states i and j. The transition rates through 
intermediate stable states can be determined by compos¬ 
ing these elementary transition rates. It is often expected 
that the most likely transition path between two stable 
states would be a direct path, but as shown below, in 
many cases it passes through intermediate stable states. 
Moreover, the transition paths and rates are generally 
asymmetric, with IF- ■ ^ R^j. In the limit of long time, 
these transitions lead to an equilibrium probability dis¬ 
tribution of occupied stable states, which we refer to as 
the limiting occupancy of the system and denote by . 
Given that we generally study large populations of cells, 
in our applications it is convenient to interpret the lim¬ 
iting occupancy as an ensemble average over different 
realizations of noise. 

To illustrate the minimum action paths and their use 
in controlling cell behavior, we first consider a two- 
dimensional model for the Caenorhabditis elegans vulval 
precursor cell (VPG) differentiation The VPG differ¬ 
entiation is representative of many other differentiation 
processes and enjoys significant experimental characteri¬ 
zation. The model describes the differentiation of VPGs 
into one of three competent lineages, 1°, 2°, and 3°, cor¬ 
responding to stable states in the system and marked as 
nodes in Fig. Ha). These lineages depend on two dimen¬ 
sionless parameters, and ^2, that determine the levels 
of epidermal growth factor (EGF) and Notch signaling, 
respectively (for the model equations, see Appendix iB]) . 
It is known experimentally that increasing EGF and de¬ 
creasing Notch signaling (relative to the base value of 
0.1) will bias cells towards lineage 1°, that decreasing 
EGF and increasing Notch will bias cells towards lineage 
2°, and that decreasing both EGF and Notch will bias 
cells towards lineage 3° [25|. Figure [2fa) shows the state 
space for the VPG system with equal, intermediate levels 
of EGF and Notch signaling = £2 = 0.1) and a real¬ 
istic level of noise {s = 0.007). The calculated transition 
paths, transition rates, and limiting occupancies are indi¬ 
cated by the edges, their width, and the size of the nodes, 
respectively. For these parameters, our calculations in¬ 
dicate that the limiting occupancy is comparable for all 
three stable states, with lineage 2° having the highest 
occupancy. 


B. Optimal Least Action Control 

We now turn to the manipulation of the tunable pa¬ 
rameters which may include, for example, gene ex¬ 
pression levels, protein activity, and interaction rates. 
Specifically, we seek to modify ft to optimize either spe¬ 
cific transition rates Rfj directly or the limiting occu¬ 
pancy through the alteration of the transition rates, 
while recognizing that the parameters can be modified 
only within specific ranges determined by biophysical and 
experimental constraints. The latter include sparsity con- 
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(a) ^1=0.1,4=0.1 (b) 4=0.24,4=0.01 
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FIG. 2. OLAC applied to a two-dimensional model of VPC differentiation, (a, b) State space representation of stable states 
(nodes) and optimal transition paths (edges) before intervention (a) and after OLAC is applied to maximize the limiting 
occupancy of lineage 1° (b). Node size indicates occupancy and edge width indicates the (negative of the) log transition 
rate along the corresponding minimum action path; the color code on the transition paths indicates the derivative of the 
quasipotential. The background shows the velocity field for the given parameters. For equal EGF and Notch signaling (a), 
lineage 2° is the one with highest occupancy. The OLAC solution (b) indicates that high EGF signaling and low Notch signaling 
will lead to the maximum occupancy of lineage 1°. (c) Trajectory in the parameter space for one realization of the optimization 
procedure, where the contour plot indicates the limiting occupancy of lineage 1°. Note that each step of the optimization 
routine increases the occupancy of this state, (d, e) NESTs for the initial (d) and optimized (e) systems. In all panels, the 
noise intensity is assumed to be £ = 0.007 (as used previously [ 13 ]). 


straints, which we can use to effectively limit the number 
of targets in the control set while avoiding a combina¬ 
torial explosion (Appendix [C]). The problem is formal¬ 
ized as the maximization of an objective function of in¬ 
terest G (unique to each problem) over the parameters 
ft subject to the given constraints {gk{^) = Ck}k and 
^ Vk^k- 

max G'({S'*/e)};£). (3) 

{hk {^)<r]k}k 

This procedure constitutes the central step of our imple¬ 
mentation of OLAC. Note that this formulation is inde¬ 
pendent of the dimension and complexity of the system. 
Thus, given a well-defined model we can identify control 
interventions able to alter the switching behavior and/or 
the stability of the states in desired ways without the 
need to know explicitly a priori how variations of the 
tunable factors affect the system (beyond the implicit 
dependence defined by the dynamical equations). 

The objective function and the associated constraints 
do not need to be linear, allowing a wide range of pos¬ 
sible dynamical behaviors to be optimized. For any set 
of such constraints, Eq. ([3]) can be solved numerically 
through a second quasi-Newton method step that nests 


the adaptive minimum action method used to determined 
the minimum action paths in connection with Eq. @. 
Importantly, OLAC is highly scalable, with the following 
computational cost: 

Cost - 0{\n\D{K (4) 

where D is the number of variables defining the state 
space, \Cl\ is the number of tunable parameters under 
consideration, P is the number of stable states, and K is 
the number of transitions upon which G depends. This 
cost estimation follows from noting that every optimiza¬ 
tion step of OLAC relies on evaluations of G, each of 
which requiring K runs of the nested optimization step, 
where each such run has a cost that is linear in D when 
using the L-BFGS optimization method [26|. Further¬ 
more, for each top-level optimization step m , every one 
of the P stable states has to be continued \fl\ times, each 
time at an integration cost that is linear in D if the av¬ 
erage network’s degree is approximately constant, as in 
most network models (it would be at most quadratic in D 
in the most general case) (Appendix [A|) . Parameter 7 ac¬ 
counts for other constants describing the relative cost be¬ 
tween optimization and integration. Because of this high 
scalability, our method can be applied to complex high¬ 
dimensional multi-parameter systems without excessive 
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computational cost. In this way OLAC expands on pre¬ 
vious foundational work that demonstrated how barriers 
between stable states can be altered to control the re¬ 
sponse to fluctuations and multistability [ 23 ^ to now 
address large networks with many variables and many 
potential control parameters. 

Before turning to high-dimensional systems, we con¬ 
sider an illustrative example application of OLAC in 
which we maximize the final occupancy of lineage 1° in 
the VPC system. As an example constraint we stipu¬ 
late that neither of the other two stable states be lost 
to bifurcation. This non-bifurcation constraint is in¬ 
tended to depict how the presence of dosing and/or ex¬ 
perimental limitations may make complete elimination 
of undesired states infeasible. This condition is im¬ 
posed as the constraint that the states = X^{ft) 
representing each lineage k remain stable fixed points: 
gk{^) = = 0 and hk{^) = Re{Xk) < —tq. In 

this notation, = Xk{^) is the eigenvalue of the Jaco¬ 
bian matrix DF{X, Q) I with largest real part, and tq is 

a tolerance (set to be 0.1 in our simulations). The objec¬ 
tive function is G{Cl] e) = vfo (Cl) and the tunable param¬ 
eters are ft = {^^ 2 }- The result, shown in Fig. Efb), 
indicates a substantial (3-fold) increase in the limiting 
occupancy of the lineage 1° state for the optimal control 
intervention identified by OLAC. The parameter-space 
path for a representative realization of the optimization 
is shown in Fig. He). The optimal intervention, defined 
as an average over multiple realizations, is a combination 
of increased EGF signaling 0.1 ^ 0.24) and reduced 
Notch signaling (^ 2 - 0.1 ^ 0.01), which has the net ef¬ 
fect of lowering the barrier for transitions from lineages 2° 
and 3° to lineage 1° while maintaining a high barrier for 
exiting this lineage. This controlled state corresponds to 
signaling strengths that have been observed experimen¬ 
tally to indeed bias cell lineages towards lineage 1° [25|. 
As mentioned above, for these results we utilized Eq. (ED 
with the proportionality constants omitted to calculate 
transition rates between states. The inclusion of these 
prefactors could potentially alter the occupancy calcula¬ 
tions of the stable states and in turn invalidate the results 
of OLAC. In the Supplemental Material [ 2 ^, Sec. SI, we 
apply sensitivity analysis to quantify the uncertainty as¬ 
sociated with omitting prefactors and show that doing so 
does not significantly alter the results in this case. That 
section also discusses more generally under which condi¬ 
tions prefactors can be omitted. 


C. Network of State Transitions 


Our formulation leads to a succinct and intuitive NEST 
representation for the transition dynamics, in which the 
nodes are unique stable states and the weighted, directed 
edges between nodes represent the rates of transition be¬ 
tween the stable states. The basis for this representation 
is the observation that, for small 5 , the trajectories of 


the system in Eq. (HD are most often close to a stable 
state or transitioning between stable states along a min¬ 
imum action path. The trajectories will rarely venture 
into other regions of the state space, allowing us to fo¬ 
cus only on this “waiting-transition” dynamics without 
sacrificing information about the system’s behavior. Eor 
a system with P stable states we can write the P x P 
transition matrix 



Ri,j i j 

i = X 


( 5 ) 


which defines a (continuous-time) Markov process on the 
stable states of the system. In our numerical calculations 
we use the fact that the limiting occupancy is described 
by the equilibrium solution of this process. In fact, the 
Markov process specifies all attributes of the associated 
NEST. The NEST is conceptually similar to the state 
transition network used in physical chemistry and bio¬ 
chemistry [ 3314 ^ as well as to transition networks stud¬ 
ied in mathematics 0. There are, however, key differ¬ 
ences between our approach and those considered previ¬ 
ously, in addition to the fact that we focus on network 
systems. In particular, the NEST does not assume the 
existence of a potential energy landscape and is defined 
for nonvanishing levels of noise, which required us to de¬ 
velop a new formulation that accounts in particular for 
long mixing times and for values of 5 larger than 5'*^ (see 
Supplemental Material [ 2 ^, Sec. S4). Eurthermore, un¬ 
like static state transition networks, the transition rates 
in the NEST are malleable and can be rationally manip¬ 
ulated with OLAC. 

Eor the VPC model, the NESTs corresponding to the 
unmodified signaling strengths [Eig. Efa)] and to the 
signaling strengths that optimize lineage 1° occupancy 
[Eig.[2jb)] are shown in Eigs.[2jd) and[2je), respectively. 
These networks represent a substantial distillation of the 
dynamics of the underlying biophysical system and can 
be used to simplify and explain the transition dynamics 
in a high-dimensional system without the need to con¬ 
sider its entire state space. In particular, for the range of 
edge widths shown, the optimized NEST in Eig.[2fe) has 
edges between all nodes except from lineage 3° to lineage 
2°, indicating that a direct transition between these two 
states is highly unlikely whereas an indirect transition is 
possible; indeed, the two-step transition 3°^1°^2° has 
an overwhelmingly higher rate (10^ times higher). By 
comparing with the NEST of the original system, shown 
in Eig. [2fd), this also demonstrates that direct transi¬ 
tions for one parameter regime can become indirect for 
another. 

The OLAC method can be implemented directly on 
the NEST representation. Indeed, the objective function 
G{{S*j}; e) is naturally defined in terms of the transition 

matrix of the system, R^. As our application to the VPC 
system shows, the combined effect of optimizing this ob¬ 
jective function is to vary the height of the transition 
barrier along the minimum action path between the sta- 
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FIG. 3. Controlled stochastic lineage switching in a multi-dimensional model of HPC differentiation, (a, b) NESTs of the 
model for the unmodified parameters (a) and for OLAC applied to optimize the limiting occupancy of the [3 cell state (red) 
(b), for £ = 0.01 In both panels, node size indicates the limiting occupancy of the state and edge width indicates the (negative 
of the) log transition rate. OLAC identifies that only three (out of ten) transcription factors need to be tuned in the model 
to optimize ^d-cell occupancy: MafA (increased), Brn4 (decreased), and (5-gene (decreased), (c) Transition hierarchy into the (3 
cell state for the optimized system. Two states ((5 and a/13) transition directly; two others (a and a/S) require passage through 
an intermediate state—in this case a/13. 


ble states. This shows that OLAC itself does not require 
determining the full quasipotential of the system, which 
would be computationally prohibitive in high dimensions. 

III. APPLICATIONS 

A. Controlling Pancreas Cell Transdifferentiation 

An important research problem in cellular reprogram¬ 
ming concerns the induction of insulin producing pan¬ 
creatic /3 cells from non-insulin producing cell lineages; 
interventions capable of achieving this goal could lead to 
new treatments for type I diabetes. In order to compu¬ 
tationally identify an optimal intervention to induce the 
desired reprogramming, we consider a ten-dimensional 
model of the hierarchical pancreas cell (HPC) differenti¬ 
ation [^. The model has five stable fixed points—three 
representing differentiated endocrine pancreas cell types 
(a, and ^) and two representing intermediate states 
{a/^ and a/5). The expression level of the ten regula¬ 
tory genes, {x^}, are each assumed to be tunable inde¬ 
pendently through a factor (details of the model are 
given in Appendix [0. 

We first consider the uncontrolled model, in which 
= 1 for i = 1,..., 10. As shown in Fig. EJa), in 
this case the two intermediate states attract the majority 
(> 99%) of the occupancy. Furthermore, transitions to 
the /3 state from the a and 6 states occur at negligibly 
small rates (< 10“^^^), indicating that such lineage re¬ 
specification effectively never occurs spontaneously. We 
apply OLAC to this model in order to identify the op¬ 
timal combination of control actions that maximize the 
occupancy of the [3 state, rj|; the admissible interven¬ 
tions are limited to those for which no bifurcations oc¬ 
cur, which is imposed using the same constraints as in the 
previous example. The optimal intervention that maxi¬ 
mizes is a three-gene one, consisting of the downregu- 
lation of Brn4 and ^-gene combined with the upregulation 
of MafA. The resulting NEST for 5 = 0.01 is shown in 


Fig. 00). Under this optimal intervention, the limiting 
occupancy of the /3 lineage goes from less than 0.01 to 
more than 0.99. 

The analysis specifically shows the reliance of some, 
but not all, lineages on two-step transitions to reach the 
desired [3 lineage [Fig. [3fc)]. Previous research has sug¬ 
gested that indirect lineage respecifications might be sub- 
optimal reprogramming strategies [38| . This is clearly the 
case for the 6^/3 reprogramming, which is optimized 
through a direct transdifferentiation event. For the a 
and a/6 lineages, however, with overwhelming likelihood 
the transformation to the /3 state will pass through the 
intermediate a/P state. Thus, which of the two cellular 
reprogramming strategies (direct or indirect) is optimal 
is context-dependent and cannot be determined without 
specification of the system, the initial state, and the fi¬ 
nal state. Systematically accounting for such context- 
dependence could lead to new advances in the develop¬ 
ment of cellular reprogramming technologies. 


B. Predicting Anti-Cancer Therapeutic Targets 

Evasion of apoptosis is one of the hallmarks of cancer 
cells [ 39 I . As such, identifying tunable factors in the cell 
death pathway that effectively eliminate proliferative (or 
abnormal survival) cell states without harming healthy 
cells could lead to new therapeutic targets. To computa¬ 
tionally identify candidate targets, we employ OLAC to 
analyze a reformulation of the Boolean model of the cell 
death pathway proposed in Ref. [40| . This reformulation 
is a continuous-variable model generated using the Hill- 
Cube methodology [4l| , which is more amenable to anal¬ 
ysis and preserves all relevant properties of the original 
model, including the stable states. The model is com¬ 
prised of 22 genes central to the programmed cell death 
and 42 parameters representing kinetic constants of the 
different reactions, which we denote by q, i = 1,.. .42 
(Appendix [Bj . It follows that two stable states form for 
all values of the parameters: apoptosis and necrosis. For 
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FIG. 4. Optimal therapeutic interventions for the cell death regulatory network model. The network has 22 genes (circles), 3 
input parameters (top), and 42 target tunable parameters (edges), where the edge heads distinguish between excitatory (arrow) 
and inhibitory (bar) regulatory interactions. OLAC is applied to induce a bifurcation transition from survival to apoptotic 
states. The parameters involved in one or more optimal intervention are consistently upregulated (green) or downregulated 
(red) by the interventions; those not recruited by any optimal intervention are shown in black. The edge width indicates the 
prevalence of that parameter in optimal interventions for the various cell types and intervention strengths considered (up to the 
elimination of the survival state). The size of each circle represents the sensitivity of the gene to pro-apoptotic interventions, 
defined as the change in gene expression between the uncontrolled scenario and the smallest-intervention scenario at which the 
survival state is eliminated; the color distinguishes up- and down-expression. Genes in white are those that change expression 
substantially between the survival and apoptotic states. These results are based on 6 simulated cell types and 9 different 
intervention strengths (Supplemental Material [2^, Fig. S3). 


a specific range of parameters representing healthy cells, 
a third state is also stable: the so-called naive state. For 
a different parameter range, however, a different stable 
state arises, corresponding to a survival cell type that 
is resistant to the apoptotic signal. This survival state 
represents cancer and is the focus of our discussion. 

Our goal is to predict therapeutic targets that can in¬ 
duce transition from the survival state to the apoptotic 
state, without increasing the rate of apoptotic death in 
normal cell types or causing them to become survival 
cells. Although noise can in principle induce switches 
from the survival to the apoptotic state, only a frac¬ 
tion of the cells would transform as desired when both 
stable states exist. We therefore neglect the effect of 
noise temporarily and show that even in this case OLAC 
can be used to identify successful optimal interventions, 
which under these conditions lead to a bifurcation that 
completely eliminates the undesired (survival) state. To 
avoid inflammatory response, it is also desirable not to 
induce transitions to the necrotic state. These condi¬ 
tions are assured by taking G{Cl) = —S'*(survival ^ 
apoptosis) as the objective function to be maximized, 
and by imposing the constraints AS* (naive ^ apoptotic) 
> —'i^o, AS* (naive ^ necrotic) > and S* (survival 
^ necrotic) > where A indicates change under the 
control intervention and is taken to be 0.05 in our sim¬ 


ulations. To encompass the largest possible set of can¬ 
didate targets, we assume that any of the 42 parameters 
of the model can be tuned in experiments, and hence we 
take ft = {ci}f^i. However, since existing experimental 
techniques cannot be easily used to manipulate a large 
number of targets, we further constrain each control in¬ 
tervention to only involve a relatively small number of all 
tunable parameters. This can be achieved by stipulating 
that the sum of the absolute values of all changes to the 
parameters must be equal to a pre-specified intervention 
strength xo, as detailed in Appendix [Cl 

We thus apply OLAC to the cell death model for the 
objective function and constraints above. To account for 
genetic heterogeneity between cancer cells, we analyze 
six different survival cell types, represented here as six 
sets of unique values for the parameters ft. These pa¬ 
rameters represent nondimensional kinetic constants for 
each gene-gene interaction. In our simulations the indi¬ 
vidual uncontrolled parameters for these cell types are on 
average 0.50 ± 0.03 and the intervention strength xo is 
taken to be 0 .1. 0.2, ... 0.9 (Fig. S3, in the Supplemen¬ 
tal Material [^, shows the breakdown of all cases). The 
number of parameters modified by optimal interventions 
tends to increase as xo increases, ranging from an aver¬ 
age of 2.7 to 5.7 for xo varied from 0.1 to 0.9. For each 
cell type, a successful intervention (eliminating the sur- 
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vival state) is always achieved for large enough xo within 
this interval. On average, approximately only 5 out of 
the 42 parameters need to be modified in the optimal 
successful interventions. To put this result in perspec¬ 
tive, we note that if the sparsity constraint in Eq. m 
is disabled, OLAC leads to an average of no less than 40 
modified parameters for a successful intervention. This 
constraint, which generalizes immediately to any system, 
is therefore effective to restrain the number of control 
parameters. 

Figure m summarizes the results, showing that a unique 
subset of only 10 parameters is needed to form the (on av¬ 
erage) 5-parameter successful target sets for any cell type. 
The biological functions and prevalence of these targets 
are explicitly indicated in Table SI of the Supplemen¬ 
tal Material [2^. Notably, only two targets are included 
in interventions found for all six cell types, namely the 
parameters whose predicted increase will decrease the ac¬ 
tivation of NF/^B by IKK and the activation of IKK by 
RIPlub (Fig.m parameters 26 and 36, respectively). The 
identification of these two targets is not entirely surpris¬ 
ing since NF/i:B is a central regulator of the cell death 
pathway whose over-activation has been implicated in the 
cellular transition to cancer [i^; the consistent identifi¬ 
cation of both targets across all six cell types is an indica¬ 
tion of the robustness of our approach. Aside from these 
global targets, OLAC also predicts unique combinations 
of targets for each cell type, in many cases indicating 
genes and interactions that have only recently been iden¬ 
tified as possible cancer targets—e.g., the potential of 
suppressing the activation of cFLIP by NF/^B [i^ (Fig. ID 
parameter 9). The identification of optimal target combi¬ 
nations that are unique to different cell types illustrates 
the potential of OLAC to assist the development of per¬ 
sonalized therapeutic strategies as well as of interventions 
to address various forms of cancer Q and to manipulate 
heterogeneous multistable cells in general 

OLAC finds an optimal control action whether or not 
a bifurcation has been reached, allowing its efficacy and 
possible adverse effects to be monitored in experiments 
as the strength of the intervention is increased. Theoret¬ 
ically, the efficacy of an intervention can be defined as 
the relative reduction of the action associated with the 
transition from the survival state to the apoptotic state 
[Supplemental Material [2^, Fig. S3(c)]. Experimentally, 
the efficacy can be more easily estimated by monitoring 
the predicted gene expression changes induced by the in¬ 
terventions [Supplemental Material [2^, Fig. S3(b)]. As 
shown in Fig. |4]for interventions that eliminate the sur¬ 
vival state, the sensitivity to the control interventions can 
vary widely across different genes. For example, in every 
cell type considered, the expressions of CASP3, SMAC, 
and CYT-c (among others) are predicted not to change 
at all until the elimination of the survival state; in con¬ 
trast, cFLIP, IKK, and NF/^B change expression in all 
cell types. 
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FIG. 5. Transition time for the cell death model as a function 
of barrier height and noise strength. Transition time was cal¬ 
culated using Eq. The black line separates those barrier 
height/noise strength combinations that occur over a thera¬ 
peutically relevant time scale (bottom right) and those that 
do not (top left). 


C. Biophysically Relevant Transition Times 

In cases where OLAC identifies an intervention to in¬ 
duce a bifurcation to eliminate an undesired state, it is 
expected from experimental work in deterministically re¬ 
programming somatic cells to a pluripotent state (a pro¬ 
totypical example of an induced cell state transition) that 
such a transition will occur on a relatively short time 
scale, within approximately one week [48|. However, due 
to constraints on parameter changes, it may not always 
be possible to induce a bifurcation. To maintain bio¬ 
logical relevance, the time scale over which these non- 
bifurcative interventions act should not be exceedingly 
long. Given a noise strength 5 and a barrier height 
between two states, the approximate mean first exit time 
flj (fl) from that state can be estimated as M 

{n) oc exp , (6) 

where the transition time increases for lower noise levels 
and higher action barriers, as expected. Since Tlj in 
Eq. (|6]) is dimensionless, this quantity is best interpreted 
as a relative increase in transition time over the case of 
Sij = 0, which from above we take to be one week. 

Figure [5] shows the mean first transition time using as a 
model application the cell death model considered in the 
previous section. The figure shows the mean transition 
time as a function of both noise level and barrier heights 
for a single cell type (as defined by each intervention 
strength in Fig. S3(c) of the Supplemental Material 13). 
It follows that for cases where 5'*^ (fl) < 3^, the average 
transition time will be less than 20 weeks. This time scale 
is a reasonable upper limit for the biological relevance of 
any intervention: because transition times are exponen¬ 
tially distributed, interventions at this strength will cause 
a measurable fraction of the population to transition in 
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just a couple of weeks. This benchmark thus provides 
a straightforward criterion to determine if an identified 
intervention will be relevant in practice. Figure [5] also 
demonstrates the important fact that it is not the size 
of the dynamical system that determines the switching 
time between states, but rather the ratio of the barrier 
height to the strength of the noise. 

We expand on this approach in the Supplemental Ma¬ 
terial Sec. S3, where we demonstrate that OLAC 
can be used to identify constrained interventions that 
achieve a desired limiting occupancy within a prespeci¬ 
fied time frame. In that section we also generalize OLAC 
to systems with (multiplicative) noise that depends on 
the system state and to systems that are best modeled 
using a chemical kinetic formulation jiiMl- 


IV. DISCUSSION 

The method proposed here, optimal least action con¬ 
trol (OLAC), represents a new direction in the control 
of biophysical systems. Traditional control approaches 
for biophysical systems are based on manipulating the 
trajectories of the system, while our method is based in¬ 
stead on manipulating the epigenetic landscape through 
which these trajectories travel. This is achieved by ef¬ 
fectively altering the barriers between different stable 
states in the quasipotential of the system, which could 
be achieved biolo gica lly through, for example, a genome 
editing approach [^. In this way OLAC stands in con¬ 
trast to other orthogonal methods that seek to control 
cellular states by directly modifying gene expression lev¬ 
els through, for example, siRNA strategies [53l-[55|. As 
such, our method naturally accounts for cellular noise 
and the incorporation of constraints on the possible con¬ 
trol actions. But in contrast to previous work that has 
sought to construct the entire quasipotential for a dy¬ 
namical system [56| , we utilize the associated network of 
state transitions (NEST) to describe and control the sys¬ 
tem at a substantially reduced computational cost, which 
renders the approach applicable to a wide range of bio¬ 
physical systems—including high-dimensional networks. 

The cell death model example illustrates one of the 
key strengths of OLAC: its effectiveness when the prob¬ 
lem involves an explosion in the number of possible re¬ 
programming combinations. In practical applications to 
multi-parameter systems, it is of interest to identify in¬ 
terventions that are not only optimal but also sparse— 
i.e., that target only a few out of many possible tunable 
factors. Such sparsity is desired since most biologically 
realizable interventions are able to control only a few pa¬ 
rameters and, for example, would not be able to directly 
control all parameters in the cell death model. Because 
OLAC can benefit from the framework of convex reg¬ 
ularization, however, this combinatorial increase in the 
number of possible interventions can be dealt with eas¬ 
ily by incorporating a sparsity constraint that has only 
marginal impact on the computational cost. 


Another key property of the approach is its flexibil¬ 
ity. For example, previous modeling approaches to ex¬ 
plain reprogramming experiments have focused m ainly 
on bifurcations that destabilize or eliminate states 
which are comparatively larger changes in the dynam¬ 
ics of the system. These approaches do not benefit from 
the presence of noise and tacitly assume that, to reach 
the desired state, the initial state must become determin¬ 
istically unstable or disappear, which, as demonstrated 
for the cell death model, may not be possible given a 
stringent enough set of constraints on the possible set of 
interventions. On the other hand, for being based on the 
Freidlin-Wentzell action, our method is effective as a uni¬ 
fied approach both (i) to exploit the presence of noise for 
control in the absence of any bifurcation (as shown for 
the pancreas model and the cell death model with low 
intervention strengths) and (ii) to identify control inter¬ 
ventions mediated by bifurcations in the absence of any 
noise (as shown for the cell death model with large inter¬ 
vention strengths). Therefore, instead of facing reduced 
performance in the presence of noise, which is a com¬ 
mon drawback of other control approaches [53|, OLAC 
benefits from the presence of noise, utilizing noise as an 
additional control tool. 

Through the use of the approach introduced here we 
have shown that, counterintuitively, the optimal lineage 
respecification trajectory is often indirect; that is, they 
correspond to cases in which the most likely trajectory 
for an optimized transition between two states will pass 
through one or more intermediate stable states. Such 
cases cannot be anticipated by common sense since for 
a therapeutic intervention, for instance, an indirect path 
may cause the cells to worsen before they improve. This 
result suggests a new possible method for identifying en¬ 
hanced reprogramming strategies, namely by systemati¬ 
cally exploring combinations of intermediate transitions. 

Our approach can also be applied to a much broader 
range of biophysical problems than those discussed in de¬ 
tail here. In particular, OLAC could be used in the con¬ 
text of synthetic biology, as researchers seek to build ever 
more complex synthetic systems and computer-aided de¬ 
sign methods play an increasingly important role [58| . In 
that context, OLAC can identify optimal parameter tun¬ 
ing to reshape the quasipotential for the rational design 
of systems with pre-specified dynamical behavior and re¬ 
sponse to noise. The same approach can also be used to 
generate insights into epigenetic diseases and the mech¬ 
anisms that give rise to them, including the possible de¬ 
pendence of their incidence rate on external versus ge¬ 
netic factors, as well as insights into potential preventive 
measures to reduce disease risks by identifying conditions 
that increase the barriers for transitions to disease states. 

Furthermore, the foundation of OLAC in the Freidlin- 
Wentzell action means that the applications of this 
method need not be biophysical. In particular, our 
method can easily be used to predict control interven¬ 
tions in other noisy multistable networks and, along with 
NEST, to characterize the basins of attraction of these 
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syste ms For example, in the Supplemental Ma¬ 

terial | 2 q|. Sec. S4, we construct the NEST for a dy¬ 
namical system with more than 100 attractors [o^ and 
demonstrate that even in a system with such substan¬ 
tial multistability, noise can effectively eliminate the oc¬ 
cupancy of the majority of attractors, leaving only a 
small fraction of them occupied. Moreover, along with 
the already mentioned applications to power-grids, poly¬ 
mer networks, and food-web networks, OLAC could also 
find use in controlling spreading processes on social net¬ 
works 0 , in inducing synchronization patterns in os¬ 
cillator networks [62|, in manipulating associative mem¬ 
ory networks [63|, and potentially in creating new at¬ 
tractors [13 . Further development of this method could 
also expand its applications to models of disease epi¬ 
demics and population dynamics. In particular, substan¬ 
tial foundational work has been done on the modeling 
of extinction events in such systems, which typically re¬ 
quires model-specific mathematical analysis j65l-l^. Be¬ 
cause a white-noise approximation, like in Eq. o ,cannot 
accurately approximate the dynamics of these systems 
(sil, expanding OLAC to apply to extinction events in 
larger networks will require using situation-specific cal¬ 
culations of the transition rates. 

Ultimately, we believe that OLAC—together with 
NEST—forms a flexible, scalable method which can be 
used to understand and control the dynamical stabil¬ 
ity and response to noise of a wide range of complex 
networks, including those with large number of dynam¬ 
ical variables, tunable parameters, and attractors. The 
method is easily implement able, with a ready-to-use nu¬ 
merical implementation included as supplemental files 
0 . The method requires no a priori information (be¬ 
yond those implicitly defined by the dynamical equa¬ 
tions) about how variations in the control parameters 
affect the system, and it can be used in concert with 
rather general constraints on the control actions. While 
OLAC can be applied to many systems, as formulated 
here application of the method requires a quantitative 
dynamical model. Extending the method to systems for 
which no model is available is an important direction 
for future research. Future work is also expected to ex¬ 
pand on the applications of the approach and to further 
demonstrate its experimental efficacy. 
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Appendix A: Calculating the Minimum Action Value 


The minimum action S^j is determined using 


Si,j = S[(pij] = mm (5 

${Ti)=ai 

${T2)=aj 


(Al) 


"1^1 = 5/. 


T 2 

Ti 




dt. 


where di and dj denote the coordinates of the initial and 
final stable states, respectively. To solve this optimiza¬ 
tion problem we minimize the discretized version of the 
functional [ 22 |, given by 




4>o,...,^ 


n=l 


AU 


-^(^n- 1 / 2 ;^) 


Atn, (A2) 


where we use Ti = to < ti < • • • < t^^ = T 2 , At^ = 

tn - tn-l, = 0(^n), and ^n-1/2 = (^n + ^n-l)/2. In 

our simulations we set T 2 = —Ti = 20 and rris = 100, 
and verified that larger values of T 2 — Ti and rris would 
not improve accuracy noticeably. 

To maximize efficiency, we regularly remesh the path 
from the time domain to the space domain and adaptively 
redefine tn according to 

rT2 

Aan / W{t)dt = W{tn-l/2)^tn, (A3) 

JTi 


where Ao^n = l/m^, tn-1/2 = (^n + tn-i)/ 2 , and 

w{t) = ^1 + C^ll^(t)|p is the monitor function mea¬ 
suring the speed along the path. We used the initial tn 
evenly spaced and C = 10^^. These calculations require 
an initial path between the stable states. The numeri¬ 
cal results we report are obtained using a straight line as 
the initial path. We have checked that using different ini¬ 
tial paths, such as those generated through the Brownian 
bridge approach [69| , the simulation always converges to 
the same optimal paths. 

We note that, since the parameters of the system are 
modified at each step of OLAC, robust numerical con¬ 
tinuation of the equilibria is necessary. We use a simple 
homotopy method to continue the stable states from ini¬ 
tial parameters U' to terminal parameters U". Specifi¬ 
cally, we use linear interpolation between these parame¬ 
ter values combined with a Newton step [^. The latter 
includes checking at each step that the desired states are 
not lost to unintended bifurcations. 

All optimization procedures are done employing the 
interior-point algorithm in the f mine on function of the 
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optimization toolbox of MATLAB . A MATLAB im¬ 
plementation of OLAC in included as supplemental files 
(see Supplemental Material [ 2 ^, Sec. S2, for details). 


Appendix B: Equations of the Models Considered 


VPC Differentiation Model. The deterministic compo¬ 
nent of the VPC model [ 2 J is given by 

dr 1 , 

Tt = 

where r = {x,yY, cti = tanh(||/+ m||)/ = 

\2x — 2 c 2 xy^ 2y -h C 2 {y‘^ — and fh = mo + iirhi + 

^ 2 ^ 2 - This is an effective representation of the combined 
effects of EGF and Notch signaling, whose signaling 
strengths are determined by and -^ 2 , respectively. The 
other parameters are r = 2.18, mo = (-0.86,-0.50)^, 
mi = (0.86, —0.50)^, and m 2 = (0.0,1.0)^. 

HPC Differentiation Model. The deterministic compo¬ 
nent of the HPC differentiation model is 
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where {as.Oe.a.kdeg, y, ym,h) = (2.2,6.0,4.0,1.0,0.25, 
0.125,4.0) are the fixed parameters in these equations. 

Each parameter cr^ represents a tunable factor to alter 
the expression of the gene represented by Xi. 

Cell Death Model. The cell death model was 
converted from the Boolean model [40| (available at 
http://www.ebi.ac.uk/biomodels-main/M0DEL0912180000) 
into a continuous version using the Odefy package 0. 

The system has 22 variables, representing gene products, 
and 42 tunable parameters. In addition, the model has 
3 input parameters that do not change in time, and 
3 output variables that indicate the state of the cell 
(distinct combinations of which indicate whether the cell 
is in the survival, apoptotic, necrotic, or naive state). 


Appendix C: Implementing the Sparsity Constraint 

In many biological systems there are dozens or hun¬ 
dreds of parameters that could potentially be changed, 
but in general only a few of them can be changed in any 
one intervention. To identify the few most promising 
targets from a large field of possible ones we employ a 
sparsity constraint. The constraint is implemented as 
|o| 

J2\An,\<xo, (Cl) 

i=l 

where, as above, A denotes the change due to the con¬ 
trol action. While the condition in Eq. (jCip is by it¬ 
self consistent with all parameters being altered, opti¬ 
mization under this constraint (as the one invoked by 
OLAC) is expected to lead to a reduced number of mod¬ 
ified parameters. The basis for this conclusion is that this 
constraint works similarly to well-established methods of 
convex regularization, which are known to lead to spar¬ 
sity under general conditions [^. The specific number 
of modified parameters as well as success rate will gener¬ 
ally depend on xo, and this dependence can be explored 
as an additional control factor. This formulation has the 
remarkable advantage of involving only one optimization 
step, and hence avoids the combinatorial explosion that 
would be involved in testing all combinations 

of possible “|fl| choose /c” tunable factors. Indeed, an ex¬ 
haustive strategy would be computationally prohibitive 
since testing for all k would require 21^1 — 1 optimiza¬ 
tion steps, which is > 10^^ for \Q\ = 42. Naturally, if 
a particular parameter is not targetable under the given 
conditions, such information can be directly be incorpo¬ 
rated into our analysis. 
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